function [counts, varargout] = getImageHist(src_img, bin_num, varargin)
% hs.img.getImageHist 直方图统计
%
% Syntax:
%   counts = hs.img.getImageHist(src_img, bin_num[ ...
%                           , 'Nodata', nan ...
%                           , 'StatisticRange', []]);
%  [counts, bin_location] = hs.img.getImageHist(___);
%  [counts, bin_location, counts_boundary] = hs.img.getImageHist(___);
%  [counts, bin_location, counts_less, counts_more] = hs.img.getImageHist(___);
%
% Params:
%   - src_img        [required ] [numeric        ] 待统计图像
%   - bin_num        [required ] [numeric;scalar ] 直方图bin个数, 记为N
%   - Nodata         [namevalue] [numeric;scalar ] 无效数据
%   - StatisticRange [namevalue] [numeric;numel=2] 统计范围, 左闭右开
%
% Return:
%   - counts Nx1向量, 表示该bin的统计个数
%   - bin_location (N+1)x1向量, 每个bin的范围(左闭右开)
%   - counts_less 标量, 表示小于统计范围最小值的像素个数
%   - counts_more 标量, 表示大于等于统计范围最大值的像素个数
%   - counts_boundary 2x1向量, 为 [counts_less, counts_more]
%
% Author: 
%   wreng, 2024年1月11日
%
% Update: 
%   wreng, 2024年1月11日
%
% Matlab Version: R2023a

%% 参数检查
arguments
    src_img (:,:) {mustBeNumeric}
    bin_num (1,1) {mustBeNumeric}
end
arguments(Repeating)
    varargin
end

narginchk(2, inf);
nargoutchk(0, 4);

%% 解析 name-value 参数
p = inputParser;
% 无效数据, 默认为 nan
addParameter(p, 'Nodata', nan, @(x)(isnumeric(x)));
% 直方图统计范围
addParameter(p, 'StatisticRange', [], @(x)(isnumeric(x)));

parse(p, varargin{:});

nodata = p.Results.Nodata;
in_range = p.Results.StatisticRange;

%% 直方图统计
% 检查bin个数
bin_num = round(bin_num);
if bin_num <= 0
    error("getImageHist errro: bin_num <= 0\n");
end

% 获取有效数据
if isnan(nodata)
    data = src_img(~isnan(src_img));
else
    data = src_img(src_img ~= nodata);
end

% bin 间隔
bin_interval = 1.0 / bin_num;

% 若统计范围为空，则默认统计所有的data
if isempty(in_range)
    in_range = [min(data), max(data)];
    in_range(2) = in_range(2) + (in_range(2) - in_range(1)) * bin_interval;
end

if in_range(1) == in_range(2)
    warning("getImageHist warn: max(img) == min(img)\n");
end

% 归一化数据
norm_data = ( data - in_range(1) ) / (in_range(2) - in_range(1));

% 转换为索引
bin_id = floor( norm_data / bin_interval );

% 小于统计范围: < in_range(1)
counts_less = 0;
% 大于等于范围: >= in_range(2)
counts_more = 0;
% 统计范围: in_range(1) <= x < in_range(2)
counts = zeros(bin_num, 1);

for i = 1:numel(bin_id)
    id = bin_id(i);
    if id < 0
        counts_less = counts_less + 1;
    elseif id >= bin_num
        counts_more = counts_more + 1;
    else
        id = id + 1;
        counts(id) = counts(id) + 1;
    end
end

% 计算每个bin的范围, 满足:  bin_location(i) <= counts(i) < bin_location(i+1)
bin_location = (0:bin_interval:1).*(in_range(2) - in_range(1)) + in_range(1);
bin_location = bin_location';

if nargout >= 2
    varargout{1} = bin_location;
end

if nargout == 3
    varargout{2} = [counts_less, counts_more];
end

if nargout >= 4
    varargout{2} = counts_less;
    varargout{3} = counts_more;
end

end
    